This notebook tests our preregistered hypothesis that there is an effect of condition on the probability of correctly recalling a studied item on a test.
library(BayesFactor)
Loading required package: coda
Loading required package: Matrix
************
Welcome to BayesFactor 0.9.12-4.2. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
Type BFManual() to open the manual.
************
library(dplyr)
Attaching package: 'dplyr'
The following objects are masked from 'package:stats':
filter, lag
The following objects are masked from 'package:base':
intersect, setdiff, setequal, union
library(data.table)
Attaching package: 'data.table'
The following objects are masked from 'package:dplyr':
between, first, last
library(forcats)
library(ggplot2)
library(purrr)
Attaching package: 'purrr'
The following object is masked from 'package:data.table':
transpose
library(readr)
library(stringr)
library(tidyr)
Attaching package: 'tidyr'
The following objects are masked from 'package:Matrix':
expand, pack, unpack
library(extrafont)
Registering fonts with R
library(wesanderson)
library(tikzDevice)
library(fst)
library(brms)
Loading required package: Rcpp
Registered S3 method overwritten by 'xts':
method from
as.zoo.xts zoo
Loading 'brms' package (version 2.12.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').
Attaching package: 'brms'
The following object is masked from 'package:stats':
ar
library(ggridges)
Attaching package: 'ggridges'
The following object is masked from 'package:ggplot2':
scale_discrete_manual
library(lme4)
Attaching package: 'lme4'
The following object is masked from 'package:brms':
ngrps
library(future)
plan(multiprocess)
# font_import() # Run once to populate the R font database with the fonts on the system
loadfonts(quiet = TRUE)
theme_poster <- theme_light(base_size = 14) +
theme(text = element_text(family = 'Merriweather Sans'),
strip.text = element_text(colour = "black"))
theme_paper <- theme_classic(base_size = 12) +
theme(axis.text = element_text(colour = "black"),
panel.grid.major.y = element_line(colour = "grey92"))
condition_colours <- wes_palette("Darjeeling1", n = 5)
condition_colours <- condition_colours[c(1,3)]
knitr::opts_chunk$set(fig.width=16, fig.height=16)
set.seed(0)
s2_test1 <- read_fst(file.path("..", "data", "processed", "s2", "s2_test1.fst"))
s2_test2 <- read_fst(file.path("..", "data", "processed", "s2", "s2_test2.fst"))
test_full <- bind_rows(s2_test1, s2_test2) %>%
mutate(block = as.factor(block),
condition = as.factor(str_to_title(condition)))
Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector
Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector
test_studied_full <- filter(test_full, studied)
Filter out excluded participant:
s2_exclude <- read_fst(file.path("..", "data", "processed", "s2", "s2_exclude.fst"))
test <- anti_join(test_full, s2_exclude, by = "subject") %>% droplevels()
Warning: Column `subject` joining factor and character vector, coercing
into character vector
test_studied <- anti_join(test_studied_full, s2_exclude, by = "subject") %>% droplevels()
Warning: Column `subject` joining factor and character vector, coercing
into character vector
According to our hypothesis, test accuracy should be higher in the fact blocks than in the default blocks, since facts are repeated with a more appropriate schedule from the start.
The plot below shows the percentage of correct answers on the test in each block, split by condition order. Only test items that appeared during the learning session are included.
test_acc_studied <- test_studied %>%
group_by(subject, block, condition) %>%
summarise(p_corr = mean(correct)) %>%
mutate(condition_order = case_when(
block == 1 && condition == "Default" ~ "Default-Fact",
block == 2 && condition == "Fact" ~ "Default-Fact",
block == 1 && condition == "Fact" ~ "Fact-Default",
block == 2 && condition == "Default" ~ "Fact-Default"
))
ggplot(test_acc_studied, aes(x = block, y = p_corr)) +
facet_grid(~ condition_order) +
geom_violin() +
geom_boxplot(width = 0.2) +
labs(x = "Block", y = "Accuracy") +
theme_poster
Make the plot shown in the paper:
p <- test_acc_studied %>%
ungroup() %>%
mutate(condition = stringr::str_to_title(condition)) %>%
ggplot(aes(y = p_corr, x = condition, group = condition)) +
geom_jitter(width = 0.1, height = 0, alpha = 0.3, aes(colour = condition)) +
geom_violin(fill = NA, width = 0.75) +
geom_boxplot(width = 0.2, outlier.shape = NA, fill = NA) +
guides(colour = FALSE) +
expand_limits(y = 0) +
scale_y_continuous(labels = scales::percent_format()) +
scale_colour_manual(values = condition_colours) +
labs(x = NULL,
y = "Test accuracy") +
theme_paper
saveRDS(p, "../output/test-accuracy.rds")
p
Make plot for presentation:
theme_presentation <- theme_light(base_size = 18) +
theme(text = element_text(family = "Merriweather Sans", colour = "black"),
strip.text = element_text(colour = "black"),
plot.title = element_text(size = 14, hjust = 0.5),
axis.title = element_text(size = 14),
rect = element_rect(fill = "transparent"),
plot.background = element_rect(fill = "transparent", color = NA),
panel.background = element_rect(fill = "white", colour = "grey70")
)
test_acc_studied %>%
ggplot(aes(y = p_corr, x = condition, group = condition)) +
geom_jitter(width = 0.1, height = 0, alpha = 0.3, aes(colour = condition)) +
geom_violin(fill = NA, width = 0.75) +
geom_boxplot(width = 0.2, outlier.shape = NA, fill = NA) +
guides(colour = FALSE) +
expand_limits(y = 0) +
scale_x_discrete(labels = c("Cold start", "Warm start")) +
scale_y_continuous(labels = scales::percent_format()) +
scale_colour_manual(values = condition_colours) +
labs(x = NULL,
y = "Test accuracy") +
theme_presentation
ggsave("../output/test-accuracy-presentation.png", device = "png", width = 4, height = 3, bg = "transparent")
As preregistered, we will evaluate the strength of evidence for an effect of condition on the probability of answering test items correctly using a logistic mixed effects model, with main effects for condition and block, and random intercepts for items and subjects. The fixed effects will be sum-to-zero contrast-coded and will have Cauchy(0,1) priors.
We will then perform a Savage-Dickey density ratio test to establish the strength of evidence for an effect of condition.
Set up contrast coding so that the model intercept reflects the grand mean p(correct), allowing us to assess the effects of condition and block separately:
contrasts(test_studied$condition) <- c(-0.5, 0.5)
contrasts(test_studied$condition)
[,1]
Default -0.5
Fact 0.5
contrasts(test_studied$block) <- c(-0.5, 0.5)
contrasts(test_studied$block)
[,1]
1 -0.5
2 0.5
Set the priors for the fixed effects:
prior_m1 <- set_prior("cauchy(0, 1)", class = "b")
Fit the model:
m1 <- brm(
correct ~ condition + block + (1 | subject) + (1 | fact_id),
family = bernoulli,
data = test_studied,
prior = prior_m1,
chains = 4,
iter = 10000,
save_all_pars = TRUE,
sample_prior = TRUE,
future = TRUE,
seed = 0,
file = "model_fits/m1"
)
Ensure that MCMC went as expected. The caterpillar plots look normal, and there is no evidence of autocorrelation in the chains.
mcmc_plot(m1, pars = c("condition", "block"), type = "trace")
No divergences to plot.
mcmc_plot(m1, pars = c("condition", "block"), type = "acf_bar")
Model summary:
summary(m1)
Family: bernoulli
Links: mu = logit
Formula: correct ~ condition + block + (1 | subject) + (1 | fact_id)
Data: test_studied (Number of observations: 2274)
Samples: 4 chains, each with iter = 10000; warmup = 5000; thin = 1;
total post-warmup samples = 20000
Group-Level Effects:
~fact_id (Number of levels: 60)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.48 0.09 0.32 0.67 1.00 8642 12770
~subject (Number of levels: 68)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.62 0.09 0.46 0.82 1.00 8443 12525
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.36 0.12 1.13 1.59 1.00 14429 14117
condition1 0.42 0.11 0.21 0.63 1.00 43086 13835
block1 0.39 0.11 0.18 0.59 1.00 45202 14866
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Visualise the model’s population-level estimates, along with their 95% credible intervals:
mcmc_plot(m1, pars = c("condition", "block"), type = "areas", prob = .95)
Warning: `expand_scale()` is deprecated; use `expansion()` instead.
The model makes the following predictions, based on its fixed effects.
# Get fitted values
new_data <- crossing(condition = levels(test_studied$condition), block = levels(test_studied$block))
fit_m1 <- data.table(fitted(m1,
newdata = new_data,
re_formula = NA,
summary = FALSE))
setnames(fit_m1, as.character(interaction(new_data$condition, new_data$block)))
fit_m1_summary <- data.table(fitted(m1,
newdata = new_data,
re_formula = NA))
setnames(fit_m1_summary, c("fit", "se", "lower", "upper"))
fit_m1_summary <- cbind(new_data, fit_m1_summary)
Overall response accuracy (across blocks and conditions):
quantile((fit_m1$Default.1 + fit_m1$Default.2 + fit_m1$Fact.1 + fit_m1$Fact.2)/4,
probs = c(.5, .025, .975))
50% 2.5% 97.5%
0.7906751 0.7524486 0.8266936
Difference in response accuracy between Fact condition and Default condition (across blocks):
fact_vs_default <- (fit_m1$Fact.1 + fit_m1$Fact.2)/2 - (fit_m1$Default.1 + fit_m1$Default.2)/2
quantile(fact_vs_default, probs = c(.5, .025, .975))
50% 2.5% 97.5%
0.06848663 0.03424775 0.10466724
Difference in response accuracy between block 2 and block 1 (across conditions):
block2_vs_block1 <- (fit_m1$Default.2 + fit_m1$Fact.2)/2 - (fit_m1$Default.1 + fit_m1$Fact.1)/2
quantile(block2_vs_block1, probs = c(.5, .025, .975))
50% 2.5% 97.5%
0.06257268 0.02829405 0.09857257
Now we can test the hypothesis that there is an effect of condition against the null hypothesis that there is not. We do this by computing the density ratio between the posterior and the prior of the condition coefficient at 0. The prior is a Cauchy(0,1) distribution, which means we can calculate its density at 0 analytically. The posterior density is estimated on the basis of the posterior samples from the model.
prior_density_at_zero <- dcauchy(0, 0, 1)
m1_posterior_samples <- posterior_samples(m1, pars = "b_condition1")$b_condition1
posterior_density_at_zero <- density_ratio(m1_posterior_samples, point = 0)
xval <- seq(-2, 2, by = 0.01)
plot_prior <- dcauchy(xval, location = 0, scale = 1)
plot_posterior <- density_ratio(x = m1_posterior_samples, y = NULL, point = xval) # Returns the density of the posterior at each point in x
m1_density <- tibble(xval = xval, Prior = plot_prior, Posterior = plot_posterior)
gather(m1_density, -xval, key = "Type", value = "density") %>%
mutate(density = ifelse(density < 0, 0, density)) %>%
ggplot(aes(x = xval, y = density, fill = Type)) +
geom_area(position = "identity", colour = "black", alpha = 0.75) +
geom_vline(xintercept = 0, lty = 3) +
annotate("point", x = c(0, 0), y = c(prior_density_at_zero, posterior_density_at_zero)) +
labs(x = "Estimate for coefficient 'condition'", y = NULL, title = "") +
scale_fill_manual(values = c("#03396c", "#b3cde0"))
ggsave(file.path("..", "output", "sd_ratio.pdf"), device = "pdf", width = 4, height = 3)
The evidence in favour of a non-zero coefficient for condition is evaluated by means of the Savage-Dickey density ratio: the ratio between the prior and posterior density at \(\beta_{condition} = 0\) (the two points on the dotted line in the plot above).
BF_condition <- prior_density_at_zero / posterior_density_at_zero
BF_condition
[1] 181.9398
The evidence in favour of an effect of condition is \(BF_{10}\) = 181.9398.
What if we choose a weaker Cauchy(0, 10) prior?
prior_m1_weak <- set_prior("cauchy(0, 10)", class = "b")
Fit the model:
m1_weak <- brm(
correct ~ condition + block + (1 | subject) + (1 | fact_id),
family = bernoulli,
data = test_studied,
prior = prior_m1_weak,
chains = 4,
iter = 10000,
save_all_pars = TRUE,
sample_prior = TRUE,
future = TRUE,
seed = 0,
file = "model_fits/m1_weak"
)
prior_density_at_zero_weak <- dcauchy(0, 0, 10)
m1_weak_posterior_samples <- posterior_samples(m1_weak, pars = "b_condition1")$b_condition1
posterior_density_at_zero_weak <- density_ratio(m1_weak_posterior_samples, point = 0)
xval <- seq(-2, 2, by = 0.01)
plot_prior_weak <- dcauchy(xval, location = 0, scale = 10)
plot_posterior_weak <- density_ratio(x = m1_weak_posterior_samples, y = NULL, point = xval) # Returns the density of the posterior at each point in x
m1_weak_density <- tibble(xval = xval, Prior = plot_prior_weak, Posterior = plot_posterior_weak)
gather(m1_weak_density, -xval, key = "Type", value = "density") %>%
mutate(density = ifelse(density < 0, 0, density)) %>%
ggplot(aes(x = xval, y = density, fill = Type)) +
geom_area(position = "identity", colour = "black", alpha = 0.75) +
geom_vline(xintercept = 0, lty = 3) +
annotate("point", x = c(0, 0), y = c(prior_density_at_zero_weak, posterior_density_at_zero_weak)) +
labs(x = "Estimate for coefficient 'condition'", y = NULL, title = "") +
scale_fill_manual(values = c("#03396c", "#b3cde0"))
Even with the much weaker prior, the evidence in favour of an effect of condition is still very strong.
BF_condition_weak <- prior_density_at_zero_weak / posterior_density_at_zero_weak
BF_condition_weak
[1] 58.95594
How about a stronger Cauchy(0, 0.5) prior?
prior_m1_strong <- set_prior("cauchy(0, 0.5)", class = "b")
Fit the model:
m1_strong <- brm(
correct ~ condition + block + (1 | subject) + (1 | fact_id),
family = bernoulli,
data = test_studied,
prior = prior_m1_strong,
chains = 4,
iter = 10000,
save_all_pars = TRUE,
sample_prior = TRUE,
future = TRUE,
seed = 0,
file = "model_fits/m1_strong"
)
prior_density_at_zero_strong <- dcauchy(0, 0, 0.5)
m1_strong_posterior_samples <- posterior_samples(m1_strong, pars = "b_condition1")$b_condition1
posterior_density_at_zero_strong <- density_ratio(m1_strong_posterior_samples, point = 0)
xval <- seq(-2, 2, by = 0.01)
plot_prior_strong <- dcauchy(xval, location = 0, scale = 0.5)
plot_posterior_strong <- density_ratio(x = m1_strong_posterior_samples, y = NULL, point = xval) # Returns the density of the posterior at each point in x
m1_strong_density <- tibble(xval = xval, Prior = plot_prior_strong, Posterior = plot_posterior_strong)
gather(m1_strong_density, -xval, key = "Type", value = "density") %>%
mutate(density = ifelse(density < 0, 0, density)) %>%
ggplot(aes(x = xval, y = density, fill = Type)) +
geom_area(position = "identity", colour = "black", alpha = 0.75) +
geom_vline(xintercept = 0, lty = 3) +
annotate("point", x = c(0, 0), y = c(prior_density_at_zero_strong, posterior_density_at_zero_strong)) +
labs(x = "Estimate for coefficient 'condition'", y = NULL, title = "") +
scale_fill_manual(values = c("#03396c", "#b3cde0"))
With the stronger prior, the evidence in favour of an effect of condition is again very strong.
BF_condition_strong <- prior_density_at_zero_strong / posterior_density_at_zero_strong
BF_condition_strong
[1] 284.4968
The analysis above did not include data from excluded participants (see Prepare data section). Verify that this does not affect the conclusion by also fitting the model to the full dataset.
contrasts(test_studied_full$condition) <- c(-0.5, 0.5)
contrasts(test_studied_full$condition)
[,1]
Default -0.5
Fact 0.5
contrasts(test_studied_full$block) <- c(-0.5, 0.5)
contrasts(test_studied_full$block)
[,1]
1 -0.5
2 0.5
m1_full <- brm(
correct ~ condition + block + (1 | subject) + (1 | fact_id),
family = bernoulli,
data = test_studied_full,
prior = prior_m1,
chains = 4,
iter = 10000,
save_all_pars = TRUE,
sample_prior = TRUE,
future = TRUE,
seed = 0,
file = "model_fits/m1_full"
)
Check the chains:
mcmc_plot(m1_full, pars = c("condition", "block"), type = "trace")
No divergences to plot.
mcmc_plot(m1_full, pars = c("condition", "block"), type = "acf_bar")
The model estimates look very similar.
summary(m1_full)
Family: bernoulli
Links: mu = logit
Formula: correct ~ condition + block + (1 | subject) + (1 | fact_id)
Data: test_studied_full (Number of observations: 2299)
Samples: 4 chains, each with iter = 10000; warmup = 5000; thin = 1;
total post-warmup samples = 20000
Group-Level Effects:
~fact_id (Number of levels: 60)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.47 0.09 0.31 0.66 1.00 7966 12855
~subject (Number of levels: 69)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.66 0.09 0.49 0.85 1.00 8900 14053
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.33 0.12 1.11 1.57 1.00 12585 14927
condition1 0.44 0.11 0.23 0.65 1.00 42349 14539
block1 0.37 0.11 0.16 0.58 1.00 43303 14300
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
m1_full_posterior_samples <- posterior_samples(m1_full, pars = "b_condition1")$b_condition1
posterior_density_at_zero_full <- density_ratio(m1_full_posterior_samples, point = 0)
plot_posterior_full <- density_ratio(x = m1_full_posterior_samples, y = NULL, point = xval) # Returns the density of the posterior at each point in x
m1_full_density <- tibble(xval = xval, Prior = plot_prior, Posterior = plot_posterior_full)
gather(m1_full_density, -xval, key = "Type", value = "density") %>%
mutate(density = ifelse(density < 0, 0, density)) %>%
ggplot(aes(x = xval, y = density, fill = Type)) +
geom_area(position = "identity", colour = "black", alpha = 0.75) +
geom_vline(xintercept = 0, lty = 3) +
annotate("point", x = c(0, 0), y = c(prior_density_at_zero, posterior_density_at_zero)) +
labs(x = "Estimate for coefficient 'condition'", y = NULL, title = "") +
scale_fill_manual(values = c("#03396c", "#b3cde0"))
The Savage-Dickey density ratio is very similar.
BF_condition_full <- prior_density_at_zero / posterior_density_at_zero_full
BF_condition_full
[1] 2068469
Fit a frequentist glmer with the same model structure. The model summary confirms that we get the same coefficient estimates as with the Bayesian model.
m1_freq <- glmer(correct ~ condition + block + (1 | subject) + (1 | fact_id),
family = binomial,
data = test_studied)
summary(m1_freq)
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: correct ~ condition + block + (1 | subject) + (1 | fact_id)
Data: test_studied
AIC BIC logLik deviance df.resid
2336.0 2364.6 -1163.0 2326.0 2269
Scaled residuals:
Min 1Q Median 3Q Max
-3.8038 0.2721 0.4266 0.5437 1.4019
Random effects:
Groups Name Variance Std.Dev.
subject (Intercept) 0.3542 0.5951
fact_id (Intercept) 0.2084 0.4565
Number of obs: 2274, groups: subject, 68; fact_id, 60
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.3474 0.1094 12.321 < 2e-16 ***
condition1 0.4273 0.1070 3.992 6.54e-05 ***
block1 0.3914 0.1065 3.675 0.000238 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) cndtn1
condition1 0.070
block1 0.016 0.014
An alternative method to test for an effect of condition is to compare the marginal likelihood of a model with the predictor to that of a model without it. We can obtain the marginal likelihoods using bridge sampling. This comparison yields a Bayes factor that can be interpreted in the same way.
Fit a second model without condition as a predictor:
m2 <- brm(
correct ~ block + (1 | subject) + (1 | fact_id),
family = bernoulli,
data = test_studied,
prior = prior_m1,
chains = 4,
iter = 10000,
save_all_pars = TRUE,
sample_prior = TRUE,
future = TRUE,
seed = 0,
file = "model_fits/m2"
)
Calculate the Bayes factor:
bayes_factor(m1, m2)
Iteration: 1
Iteration: 2
Iteration: 3
Iteration: 4
Iteration: 5
Iteration: 1
Iteration: 2
Iteration: 3
Iteration: 4
Iteration: 5
Iteration: 6
Estimated Bayes factor in favor of m1 over m2: 212.21710
Again, it’s basically the same as before.
How does \(BF_{10}\) develop as more participants are added? The following code fits the model repeatedly, each time adding several more participants, and calculates the \(BF_{10}\) for each fit.
subjects <- unique(test_studied$subject)
set_sizes <- seq(1, length(subjects), by = 1)
evidence <- tibble(n = 0, bf = NA, posterior = list(rcauchy(n = 20000, location = 0, scale = 1)))
for (i in set_sizes) {
# Subset the data
test_studied_subset <- filter(test_studied, subject %in% subjects[1:i])
# Fit the model
m1_seq <- brm(correct ~ condition + block + (1 | subject) + (1 | fact_id),
family = bernoulli,
data = test_studied_subset,
prior = prior_m1,
chains = 4,
iter = 10000,
save_all_pars = TRUE,
sample_prior = TRUE,
future = TRUE,
refresh = 0,
open_progress = FALSE,
seed = 0,
file = paste0("model_fits/m1_seq", i))
m1_seq_posterior_samples <- posterior_samples(m1_seq, pars = "b_condition1")$b_condition1
posterior_density_at_zero_seq <- density_ratio(m1_seq_posterior_samples, point = 0)
BF_condition_seq <- prior_density_at_zero / posterior_density_at_zero_seq
evidence_seq <- tibble(n = i, bf = BF_condition_seq, posterior = list(m1_seq_posterior_samples))
evidence <- bind_rows(evidence, evidence_seq)
}
The plot below shows the sequential Bayes factors resulting from this analysis. (Bayes factors calculated from very small samples sizes should obviously be treated with caution.) The accumulation of evidence in favour of an effect of condition is clear. As the BF becomes larger, there are some large upward spikes; this is because by this point there is essentially no posterior mass left at zero, which makes the SD density ratio very sensitive to small fluctuations in the estimated posterior due to the randomness of MCMC sampling.
The days of data collection are marked in the plot. The stopping rule, which was evaluated at the end of each day, required there to be data from at least 40 participants and a BF of 10 or greater in either direction (this boundary is indicated by horizontal lines). At the end of day 2 we had collected data from more than 40 participants, but the BF was still within the boundaries, so we collected data for another day.
seq_plot <- ggplot(filter(evidence, n > 0), aes(n, bf)) +
annotate("rect", xmin = 0, xmax = 27, ymin = 1/100, ymax = 1000000, fill = "#1b9e77", alpha = 0.2) +
annotate("rect", xmin = 27, xmax = 42, ymin = 1/100, ymax = 1000000, fill = "#d95f02", alpha = 0.2) +
annotate("rect", xmin = 42, xmax = 68, ymin = 1/100, ymax = 1000000, fill = "#7570b3", alpha = 0.2) +
annotate("text", label = c("Day 1", "Day 2", "Day 3"), x = c(13.5, 34.5, 55), y = 0.02) +
geom_hline(yintercept = c(0.1, 1, 10), lty = c(1, 3, 1)) +
geom_vline(xintercept = c(0, 27, 42, 68), lty = 2) +
geom_line() +
geom_point() +
scale_y_log10(limits = c(1/100, 1000000), breaks = c(1/100, 1/10, 1, 10, 100, 1000, 10000, 100000, 1000000), labels = function(n) {format(n, scientific = FALSE, drop0trailing = TRUE)}) +
labs(x = "Sample size", y = "BF10", title = "Sequential Bayes factors in favour of an effect of condition")
seq_plot
The plot below shows how the posterior for the condition coefficient develops as the sample size increases. The dashed vertical line represents the point at which the posterior density is compared to the prior density (shown at the bottom in grey).
unnest_legacy(evidence) %>%
ggplot(aes(x = posterior, y = n, group = n, fill = n == 0)) +
geom_vline(xintercept = 0, lty = 2) +
geom_density_ridges2(scale = 10, size = 0.25) +
scale_x_continuous(limits = c(-2, 2)) +
labs(x = "Estimate for coefficient 'condition'", y = "Sample size", title = "Sequential posteriors for 'condition'") +
scale_fill_manual(values = c("#b3cde0", "#636363")) +
guides(fill = FALSE)
Picking joint bandwidth of 0.027
Warning: Removed 16332 rows containing non-finite values
(stat_density_ridges).
sessionInfo()
R version 3.6.3 (2020-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.5 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=nl_NL.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=nl_NL.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=nl_NL.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=nl_NL.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] future_1.13.0 lme4_1.1-21 ggridges_0.5.1
[4] brms_2.12.0 Rcpp_1.0.2 fst_0.9.0
[7] tikzDevice_0.12 wesanderson_0.3.6 extrafont_0.17
[10] tidyr_1.0.0 stringr_1.4.0 readr_1.3.1
[13] purrr_0.3.2 ggplot2_3.3.2 forcats_0.4.0
[16] data.table_1.12.2 dplyr_0.8.3 BayesFactor_0.9.12-4.2
[19] Matrix_1.2-18 coda_0.19-2
loaded via a namespace (and not attached):
[1] minqa_1.2.4 colorspace_1.4-1 ellipsis_0.3.0
[4] rsconnect_0.8.15 markdown_1.0 base64enc_0.1-3
[7] listenv_0.7.0 rstan_2.19.2 MatrixModels_0.4-1
[10] DT_0.7 mvtnorm_1.1-1 bridgesampling_0.7-2
[13] codetools_0.2-16 splines_3.6.3 knitr_1.23
[16] shinythemes_1.1.2 bayesplot_1.7.0 jsonlite_1.6
[19] nloptr_1.2.1 Rttf2pt1_1.3.7 shiny_1.3.2
[22] compiler_3.6.3 backports_1.1.4 assertthat_0.2.1
[25] cli_1.1.0 later_0.8.0 htmltools_0.3.6
[28] prettyunits_1.0.2 tools_3.6.3 igraph_1.2.4.1
[31] gtable_0.3.0 glue_1.3.1 reshape2_1.4.3
[34] vctrs_0.2.2 filehash_2.4-2 nlme_3.1-149
[37] extrafontdb_1.0 crosstalk_1.0.0 xfun_0.7
[40] globals_0.12.4 ps_1.3.0 mime_0.7
[43] miniUI_0.1.1.1 lifecycle_0.1.0 gtools_3.8.1
[46] MASS_7.3-51.4 zoo_1.8-6 scales_1.0.0
[49] colourpicker_1.0 hms_0.4.2 promises_1.0.1
[52] Brobdingnag_1.2-6 parallel_3.6.3 inline_0.3.15
[55] shinystan_2.5.0 mvnfast_0.2.5 yaml_2.2.0
[58] pbapply_1.4-0 gridExtra_2.3 loo_2.3.1
[61] StanHeaders_2.19.0 stringi_1.4.3 dygraphs_1.1.1.6
[64] boot_1.3-25 pkgbuild_1.0.3 rlang_0.4.4
[67] pkgconfig_2.0.2 matrixStats_0.55.0 evaluate_0.14
[70] lattice_0.20-41 labeling_0.3 rstantools_1.5.1
[73] htmlwidgets_1.3 tidyselect_0.2.5 processx_3.3.1
[76] plyr_1.8.4 magrittr_1.5 R6_2.4.0
[79] pillar_1.4.2 withr_2.1.2 xts_0.11-2
[82] abind_1.4-5 tibble_2.1.3 crayon_1.3.4
[85] rmarkdown_1.13 grid_3.6.3 callr_3.2.0
[88] threejs_0.3.1 digest_0.6.19 xtable_1.8-4
[91] httpuv_1.5.1 stats4_3.6.3 munsell_0.5.0
[94] shinyjs_1.0